// This header file should be included in all classes.
//
// CLASSES ==========================================
// For dynamic system simulation
// - blocks
// - simulation
// 
// For modeling
// - surface2D
// - eigenfun_2D
// - surface1D
// - eigenfun_1D
//
// for KMC simulation
// - MCArray
//
// UNIT ============================================
// Length:{nm}|layer|site
// Time: {second}
// 
// CONVENTION TO WRITE A BLOCK =====================
// == Category of variables ==
// According to the order they should appear in the declaration
// - INPUT
// - OUTPUT
// - STATE
// - PARAMETERS
// - NUMERICAL SOLUTION PARAMETERS
// - Internal auxilary variables
//
// ==============================================================
// ==============================================================


#ifndef SIGNAL_H
#define SIGNAL_H

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#include <ctime>
#include <cfloat>

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdexcept>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

using namespace std;

// CONSTANT USED BY ALL SIMULATIONS
const int    COLUMNWIDTH = 15;
const double pi = 3.141592653;
const double LPCVD_EPSILON = 1e-9; // Error control
const double SQRT_TWO = 1.414213562373095;


// --------------------------------------------------------------
/*class signal
{
int     dimension;
char    name[1024];
public:
signal(const char* iname,double it_end,double idt,int idimension,int ilog_mode=0);
~signal();
void set(double ivalue);
void set(double *ivalue);
void set(double itime,double ivalue);
void set(double itime,double *ivalue);
double get(double itime);
void reset();
};*/

// --------------------------------------------------------------
// BLOCK models a discrete time dynamic system, which has the general form
//   x_{k+1} = Ax_{k}+Bu_{k}
//   y_{k} = C_{k}x_{k}+Cu_{k}
// or more general
//   x_{k+1} = f(x_{k},u_{k})
//   y_{k} = g(x_{k},u_{k})
// --------------------------------------------------------------
class block{
protected:
  char   name[1024];
  // If time = t, then the output is valid between [time-dt,time).
  double time;      
  double dt;

  int    input_num;
  int*   input_dim;
  double** inputs;

  int    output_num;
  int*   output_dim;
  double** outputs;
  //double** outputs_next;

  // For systems without dynamics (e.g., controller, 
  // constant multipler), there is no state and the output y_{k} is 
  // a function of input u_{k}.
  //int     state_num;  
  //double* states;
  //double* states_next;
public:
  block();
  block(int i_num, const int* i_dim,
	int o_num, const int* o_dim,
	double idt);

  virtual ~block();

  bool    initialize(int i_num, const int* i_dim,
		     int o_num, const int *o_dim,
		     double idt);
  bool    ConnectInputTo(int input_ID,block& source,int output_ID);
  
  int     getOutputNum(){return output_num;};
  int     getOutputDim(int output_ID) {return output_dim[output_ID];} ;
  double* getOutputAddress(int output_ID){return outputs[output_ID];};
  int     getInputNum(){return input_num;};
  double  getDt(){return dt;};

  // Virtual functions =======================================
  
  // void reset()
  // Is called at the beginning of each single simulation
  // NOTE: 
  //   * Don't reset input !!!
  //   * Only reset output and state to their default initial value
  virtual void reset() = 0;
  
  // void update(double itime)
  // Use the state and input at itime to calculate the state at itime+dt
  // and calculate the output function
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // IMPORTANT: t_{k} = itime, t_{k+1} = itime+dt
  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  // Part I: Calculate output
  // if t_{k+2} <= itime
  //    report error
  // if itime < t_{k}
  //    report error
  // while(t_{k+1} <= itime < t_{k+2})
  //    x = x_{k+1}
  //    time = t_{k+1}
  //    Update state IF NECESSARY (see part II)
  //    x_{k+1} = Ax_{k}+Bu_{k} or f(x_{k},u_{k})
  // if t_{k} <= itime < t_{k+1}
  //    y(itime) = y(t_{k}) = Cx_{k}+Du_{k}
  //                        = g(x_{k},u_{k})
  // 
  // Part II: update state
  // 1. For discrete system with dynamics:
  //   x_{k+1} = Ax_{k}+Bu_{k}
  //   x_{k+1} = f(x_{k},u_{k})
  //   Calculate x_{k+1}
  //
  // 2. For continuous system:
  //   dx/dt = Ax+Bu
  //   dx/dt = f(x,u)
  //   Calculate x_{k+1} by integrating from t_{k} to t_{k+1}.
  //
  // 3. For system with no dynamics (states)
  //   Do nothing
  //

  // For dynamics systems (the input-output is described by an difference equation or differencial equation,
  // Update will NOT influence the output at itime.
  // For other systems, for example, P controller, constant, update will influence the output.
  // 
  // Part II: output function
  // 
  // ============== Example code: system with dynamics ================
  // void update(double itime){
  //    if(time > itime){
  //       cerr << "ERROR MESSAGE";
  //       exit(0);
  //    }
  //    while(time+dt <= itime || time == 0){
  //       The current state is out-of-date, need to update
  //       if(time != 0){
  //         x(time) = x(time+dt); // This has been calculated when 
  //       }
  //
  //       1) Read the input and update the model parameter.
  //       2) Calculate the states at time+dt;
  //       
  //       if(time==0)
  //         break;
  //       else
  //         time = time+dt;
  //	}
  //    Calculate the output. After the while loop, 
  //      time <= itime < time+dt
  //    And the output is valid for [time,time+dt). 
  //
  //    return;
  // }
  // ============= Example code: system without dynamics ==============
  // void update(double itime){
  //   if(time> itime){
  //     cerr << "ERROR MESSAGE";
  //     exit(0);
  //   }
  //
  //   while(time+dt<=itime){
  //     time = time+dt;
  //   }
  // 
  //   Calculate the output
  //   return;
  // }
  virtual void update(double itime) = 0; 
  
  // virtual void output(double itime) = 0; // Copy outputs_next[][] to outputs[][]

  // virtual set(int para_ID,double para_val);
private:
  block(const block & blk);
  block& operator=(const block& blk);
};
// --------------------------------------------------------------
class simulation{
  double dt;
  double t_end;
  time_t time_begin;
  time_t time_end;
  block  *Blocks;
  int    repeat;
  int    JobID;
  int    Mode;
public:
  simulation();
  simulation(double dt,double t_end,int iRepeat=1,int iJobID=1,int iMode = 0);
  ~simulation();
  int run();
  //void connect(block Source, int O_ID, block Destination, int I_ID);
};
// --------------------------------------------------------------
class MCSimulator{
private:
  int		EventNum;
  int*	Events;
  double* Prob;
  double  TotalProb;
  gsl_rng* rng;
  double   randnum;
public:
  MCSimulator(int iEventNum);
  ~MCSimulator();
  int TryOnce();
  void UpdateProb(int EventID,double NewProb){
    TotalProb = TotalProb-Prob[EventID]+NewProb;
    Prob[EventID]=NewProb;
  };
  void AddEvent(int EventID, double NewProb);
  double getTotalProb(){return TotalProb;};
};
// --------------------------------------------------------------
class Surface1D{
private:

public:
  int	LatticeSize;
  double* x;
  double*	h;
  Surface1D(int LatticeSize,double dx);
  ~Surface1D();
  double getR2();
  double getM2();
  double getH();
  void   getXCov(double *x);
  double getAutoCovLength();
  double getPSD();
  void   getAlphas(double *alphas);
};
// --------------------------------------------------------------
class Spectrum1D{
private:
  int     ModeNum;
  int*    ModeID;
  double* Frequency;
  double* ModeVal;	
  string  EstMethod;
public:
  Spectrum1D(int ModeNum);
  ~Spectrum1D();
  void psd(Surface1D* surface);
};
// --------------------------------------------------------------
void xcov(double* i_x, double* o_xcov_val,int L);
double ranO(int *idum);
double RandRand();
double RandN1();
unsigned time_seed();
// --------------------------------------------------------------
class MCArray{
private:
  // int array_size;
  //  int array_index_begin;
  int array_index_end;
  int size;
  int layer;
  double* array; // Probabilities
  int array_size;
  int array_index_begin;
  bool NORMALIZE_PROBABILITY;

  bool EXTERNAL_RNG; // TRUE if external random number generator is used; otherwise, each MCArray will use its own random generator
  const	gsl_rng_type *gsl_T;
  gsl_rng *gsl_r;

  inline int left_child(int parent_index){
    return 2*parent_index+1;
  };

  inline int right_child(int parent_index){
    return 2*parent_index+2;
  };

  inline int parent(int child_index){
    return (child_index-1)/2;
  }
public:
  MCArray(int isize,gsl_rng *random_num_generator = NULL);
  MCArray(const MCArray & ma);
  MCArray & operator=(const MCArray & ma);
  ~MCArray();
  //double operator[](int i);
  void   SetElement(int index,double value); // Set the probability of event #INDEX to be VALUE
  double GetElement(int index); // Return the probability of event #INDEX
  int    FindElement(double x); // Return event index
  int	 PickAnEvent();
  double GetTotal(){return array[0];};// Return the sum of all the element;

};

// --------------------------------------------------------------
class MathFunction{
private:
  int num_input;
  int num_para;
public:
  MathFunction(int i_num_input=1,int num_para = 0);
  virtual ~MathFunction();
  virtual double value(const double *x)=0;
  virtual void gradient(const double *x,double *grad)=0;
  virtual bool parameterize(const double *para);
  virtual bool parameterize(int para_ID,double para_val);
  virtual bool parameterize(int para_ID,double *para_val);
  int get_num_input(){return num_input;};
  void plot(const double *x_min,const double *x_max,const double *dx,string FileName);
};

// Utilities function to test (derived class) of  MathFunction
int perturb_MathFunction(MathFunction* fun,const double *x,const double *epsilon);
void MathFunction_plot(MathFunction* fun, 
		      const double *x_min,
		      const double *x_max,
		      const double *dx,
		      const int dimension,
		      bool (*FuncPtrValidRegime)(double* x,int dimension),
		      const double VALUE_FOR_INVALID_REGIME,
		      string FileName);
// --------------------------------------------------------------
#endif

